.. _tutorial_7:
=======================================================
TUTORIAL 7 : Introduction to latticeswitch Monte Carlo
=======================================================
:Author:
Tom L. Underwood, t.l.underwood{at}bath.ac.uk
Introduction
============
A common problem is as follows: given two solid phases, which is more stable at a given temperature and pressure? The more stable phase is
that with the lower free energy of the two. Thus calculating the *free energy difference* between the two phases reveals which phase is more stable.
Latticeswitch Monte Carlo (LSMC) is a powerful method for calculating such free energy differences. Here we demonstrate
how to use LSMC in DL_MONTE to calculate the free energy difference between two solid phases at a given temperature and pressure.
Specifically, we use LSMC to investigate which of the fcc or hcp crystal phases of the LennardJones system are stable on the :math:`P=0` isobar.
The question of whether the fcc or hcp phase is the stable phase of the LennardJones solid has a long history, and was resolved only relatively
recently. The answer is that the hcp phase is stable at lower temperatures and pressures, while the fcc phase is stable at higher temperatures and
pressures (below the melting line).
We will use LSMC to investigate whether this behaviour is borne out in simulations using 216particle
cell and the following incarnation of the LennardJones system: each particle in the cell interacts with all other particles via the
LennardJones potential
.. math::
\phi(r) = 4\epsilon\biggl[\Bigl(\frac{\sigma}{r}\Bigr)^{12}\Bigl(\frac{\sigma}{r}\Bigr)^{6}\biggr]
with no cutoff, using the minimum image convention (and there is no other contribution to the total energy, e.g. due to tail corrections).
Prerequisites

This tutorial assumes that the reader is familiar with the basics of DL_MONTE, i.e. can perform simple NVT and NPT simulations.
Moreover it is also assumed that the reader has a basic understanding of biased sampling techniques.
Solutions

Note that example output and input files for all simulations to be performed in this tutorial can be found in the *solutions* directory.
This is important in light of the fact that some of the simulations could take a long time to run: one may wish to skip the actual running
of the simulation and proceed straight to the analysis or inspection of the output files.
Background  LSMC
=================
We begin by providing a short description of LSMC. More detailed information can be found in the user manual, as well as the references provided at
the end of the tutorial.
Consider the free energy difference between two phases 1 and 2. This is given by the following equation:
.. math::
\Delta G \equiv G_1  G_2 = k_BT\log(p_1/p_2)
where :math:`p_1` is the probability being in phase 1 and :math:`p_2` is the probability of the system being in phase 2 (assuming that the system can
only be in either phase 1 or phase 2); :math:`G_1` is the free energy of phase 1, and similarly for :math:`G_2`;
and :math:`T` is the temperature of the system. Note that we are treating a system at constant temperature and constant pressure here, so the free energy
here is the *Gibbs* free energy. In principle the above equation could be exploited to calculate :math:`\Delta G` from simulation as follows:
measure the times :math:`t_1` and :math:`t_2` that the system spends in phases 1 and 2 during the simulation, and then use the fact that
:math:`t_1/t_2=p_1/p_2` in the above equation to obtain :math:`\Delta G`.
Of course this approach requires that the both phases are explored on timescales accesible to simulation. Unfortunately this is not the case in
conventional methods which invoke 'realistic' particle dynamics. In such methods the time taken for the system to transition from one phase
to the other and back again is simply too long for the above approach to be tractable.
The key aspect of LSMC is a *lattice switch* Monte Carlo move which, if accepted, takes the system *directly* from the current phase to the 'other'
phase. In theory, if lattice switch moves are attempted and accepted frequently, then the system would explore both phases on reasonable timescale.
This in turn would enable :math:`\Delta G` to be evaluated via the above equation.
In a lattice switch move a configuration corresponding to the other phase is constructed from the current configuration. The new configuration is
constructed by 'switching' the underlying lattice which characterises the current phase for a lattice which characterises the other phase,
while preserving the *displacements* of all particles from their lattice sites. This is best illustrated pictorially; see the figure below.
.. figure:: ImagesPSMC/tutorialpsmc_latticeswitch.jpg
:width: 600px
Illustration of a lattice switch move. The current configuration (left) corresponds to the square phase, and the move will take the system
to a new configuration in the triangular phase (right). Note that in the current configuration the particles are all close to their lattice sites
(red crosses, which form a square lattice). The move transforms the underlying square lattice (red crosses) into a triangular lattice (blue
crosees), while keeping the displacements of the particles from their lattice sites (black arrows) unchanged  the displacement of particle
:math:`n` is the same in the current and new configurations.
However this is not the whole story. While lattice switch moves provide a means for generating reasonablelooking configurations in the other phase, these
configurations are usually of very high energy, and accordingly the Metropolis algorithm almost always rejects lattice switch moves. Thus lattice
switch moves alone are not enough to enable both phases to be explored in a single simulation. The extra ingredient which is required is *biased
sampling*, which we exploit to increase the proportion of lattice switch moves which are accepted to an acceptable level. To elaborate, we first define
an order parameter :math:`M` which distinguishes the following three types of configuration:
1. configurations realised in phase 1 at equilibrium; those which one would see in a conventional Monte Carlo simulation of phase 1. These configurations
have large negative values of :math:`M`, i.e. :math:`M\ll 0`. We refer to these configurations as phase1 *equilibrium configurations*
2. configurations realised in phase 2 at equilibrium; those which one would see in a conventional Monte Carlo simulation of phase 2. These configurations
have large positive values of :math:`M`, i.e. :math:`M\gg 0`. We refer to these configurations as phase2 *equilibrium configurations*
3. configurations where lattice switch moves *both from, and to* the configuration are highly likely. As described in a moment, these configurations act
as a gateway between the two phases, for which reason we refer to them as *gateway configurations*. They have small values of :math:`M`,
i.e. :math:`M\approx 0`
By defining :math:`M` in this way we have created a pathway linking the two phases; traversing :math:`M` takes us between phases.
Consider a phase1 equilibrium configuration. Here, :math:`M\ll 0`, and lattice switch moves are fruitless  they are essentially always rejected.
Gradually increasing :math:`M` to :math:`\approx 0` however takes the system to gateway configurations where lattice switch moves are likely to be
accepted, and the system can transition (via a lattice switch move) to phase 2. Increasing :math:`M` further
eventually takes us to the phase2 equilibrium configurations, where lattice switch moves are again fruitless.
With this in mind we perform a simulation using biasing to explore the whole range of :math:`M` uniformly. The result is that both the equilibrium
configurations for each phase, and the gateway configurations which are necessary for sucessful switching to and from both phases, are sampled.
There is of course the question of how :math:`M` is defined so that it has the desirable properties just mentioned. This is beyond the scope of
the current tutorial; it is not esssential to know this to effectively use LSMC in DL_MONTE. However the references
provided at the end of the tutorial for further information regarding how :math:`M` is defined.
LSMC workflow

We now turn to the practicalities of LSMC free energy calculations. Calculating a free energy difference using LSMC in fact involves a number of steps:
1. A *preliminary simulation* for each phase.
Above we said that we would use biased sampling to sample the whole range of :math:`M` uniformly. Of course, in practice we must choose upper and lower
limits to the range of :math:`M` which our simulation will explore. However, we do not know *a priori* what limits are appropriate. We must
therefore perform a short preliminary simulation for each phase to determine the appropriate range of :math:`M` to use in subsequent simulations.
These simulations should be 'conventional' Monte Carlo simulations, i.e. lattice switch moves and biasing should not be invoked. However the
LSMC order parameter :math:`M` must be tracked during the simulation.
2. A *biasgeneration simulation*.
Using biased sampling to sample uniformly over the range of :math:`M` (determined in the preliminary simulations) requires knowledge of the bias
function which results yields uniform sampling. Unfortunately this too is not known *a priori*. Accordingly the next
step is to perform a biasgeneration simulation to *learn* the bias function which leads to uniform sampling over our chosen range of :math:`M`.
3. A *production simulation*: a biased simulation using the bias function determined in the last step, which samples uniformly
over :math:`M` in the chosen range, and accordingly explores both phases. From the data obtained from this simulation we can determine
:math:`p_1/p_2`, and hence :math:`\Delta G` using the equation given at the beginning of this section.
Below we provide a walkthrough of this procedure for calculating the hcpfcc free energy difference of the LennardJones system described at the
beginning of this tutorial at temperature :math:`T=0.1\epsilon/k_B` (where :math:`\epsilon` is the LennardJones well depth and :math:`k_B` is
Boltzmann's constant) and pressure :math:`P=0`, using DL_MONTE.
Preliminary simulations
=======================
A basic LSMC simulation in DL_MONTE requires 5 input files: CONTROL, FIELD, CONFIG, CONFIG.1 and CONFIG.2.
The first 3 of these files are required in all DL_MONTE simulations. However CONFIG.1 and CONFIG.2 are specific to
LSMC. What follows is a description of the input files for our first preliminary simulation. These files can be found in the
main tutorial directory.
FIELD

We begin with the FIELD file (note that the line numbers at the start of each line are not actually present in the file):
.. codeblock:: html
:linenos:
LennardJones potential with sigma=epsilon=1, cutoff infinite, and 2 configs for PSMC, energy units in k_BK
CUTOFF 1000.0
UNITS k
NCONFIGS 2
ATOMS 1
LJ core 1.0 0.0
MOLTYPES 1
lj
MAXATOM 216
FINISH
VDW 1
LJ core LJ core lj 1.0 1.0
CLOSE
The FIELD file one would use in a LSMC simulation is almost the same as one would use in any other simulation. The
exception is that *NCONFIGS* must be set to 2 (see line 2 above). This is because two configurations are always used by DL_MONTE in a LSMC
simulation. The reasons for this are technical, and explained in the user manual.
There are two somewhat esoteric aspects to the FIELD file above worth remarking upon. Firstly,
we have set *CUTOFF* to 1000.0. This amounts to applying no cutoff to the LennardJones potential  as we said we would do at the beginning
of this tutorial.
Secondly, we are using an energy unit of :math:`k_B` K (where K is the kelvin). With this in mind, the LennardJones parameter
:math:`\epsilon` (the first '1.0' on line 12) is set to 1 :math:`k_B` K. Moreover the parameter :math:`\sigma` (the second '1.0' on line 12)
is set to 1 angstrom. This choice of units is for convenience. Recall that we wish to consider a temperature :math:`T=0.1\epsilon/k_B`. This
temperature, given that :math:`\epsilon=1k_B` K, corresponds to :math:`T=0.1` K  which is reflected in the CONTROL file provided in a moment.
CONFIG, CONFIG.1 and CONFIG.2

We now consider the configuration files CONFIG, CONFIG.1 and CONFIG.2. As alluded to above, CONFIG is expected to contain two
configurations in LSMC. These are prospective starting configurations for each phase. If the simulation is to be started in phase 1
(via a directive *initactive* in the CONTROL file  described in a moment), then the first configuration in CONFIG will be
used as the starting configuration for the simulation. On the other hand if the simulation is to be started in phase 2, then the second configuration
in CONFIG will be used as the starting configuration. By contrast CONFIG.1 and CONFIG.2 each contain only one configuration. The purpose
of CONFIG.1 and CONFIG.2 is to define the underlying lattices to be used in lattice switch moves: the cell and particle positions in
CONFIG.1 define the lattice for phase 1; and the cell and particle positions in CONFIG.2 define the lattice for phase 2. Details of
exactly how the lattice switch move is constructed from the information in CONFIG.1 and CONFIG.2 can be found in the user manual.
Here we choose phase 1 to be hcp and phase 2 to be fcc. Thus CONFIG.1 and CONFIG.2 correspond to perfect (i.e. undistorted) hcp and fcc crystals
respectively. Moreover we choose to start the simulation with the particle positions forming a perfect crystal. Thus the first configuration
in CONFIG corresponds to a perfect hcp crystal, and the second configuration in CONFIG corresponds to a perfect fcc crystal. Of course
one is free to use 'equilibrated' configurations, where the particle positions form an approximate crystal lattice, as prospective starting configurations
in CONFIG.
CONTROL

Below is the CONTROL file. *Note that in DL_MONTE latticeswitch Monte Carlo is referred to, for reasons we won't discuss here, as phaseswitch
Monte Carlo (PSMC). You should regard the terms PSMC and LSMC as synonymous here.*
.. codeblock:: html
:linenos:
PSMC simulation LJ # Comment line
use fed psmc # Use PSMC in 'FED' functionality, followed by block of PSMC flags; begins FED block
switchfreq 0 # Latticeswitch move frequency; NO SWITCH MOVES
initactive 1 # Phase to start in (1 or 2; chooses configuration in CONFIG)
psmc done # End of PSMC block
fed method ee 1.0 1.0 1000000000 # Use EE method for learning the bias function; NEVER UPDATE DURING SIM; NO BIASING!
fed order psmc 100 1.0E10 1.0E10 1 # Use PSMC order parameter for biasing; USE INFINITE RANGE
fed done # End of FED block
finish # End of 'use' block
ranseed # Use a random seed
steps 320000 # Simulation length
temperature 0.1
pressure 0.0
sample coordinate 21600 # Frequency to store configurations during simulation
archiveformat dlpoly4 # Configurations during simulation stored in DL_POLY4 format (HISTORY)
yamldata 216 # Frequency to store energies etc. in yaml format (YAMLDATA)
move atoms 1 216 # Move atoms, followed by frequency
LJ core
move volume ortho log 1 # Move volume, followed by freqency
start
The important features of this file are as follows. Firsty, we have used the *use fed psmc* directive (line 3), which brings
the LSMC machinery within DL_MONTE into play. This directive doubles as the start of the 'PSMC block' for specifying directives specific
to LSMC, the end of which is signified by the *psmc done* directive (line 8). Here we see that there are two directives in this
block, *switchfreq* and *initactive*. These are the only compulsory directives within the PSMC block; there are other optional
directives for controlling various aspects of LSMC within DL_MONTE, but these are not required here. *switchfreq* specifies the
frequency with which lattice switch moves are to be attempted. We have set this to 0 here because lattice switch moves are not
required yet: recall that the goal of our preliminary simulations is to determine the range of order parameter exhibited by
the equilibrium configurations of each phase, for which lattice switch moves are unnecessary (and, in fact, counterproductive:
enabling lattice switch moves at this point would result in the system immediately jumping into the phase of lowest energy and
remaining there during the preliminary simulation, as opposed to exploring the particular phase we are interested in). The
other compulsory directive is *initactive*, which specifies which phase to start in; the corresponding configuration in CONFIG
is used as the starting configuration.
Moving beyond the PSMC block, we come to the two compulsory FED directives, *fed method* and *fed order [param]* (lines 10
and 11). Note that in the latter we have specified, via *fed order psmc*, that the order parameter be the 'PSMC' order parameter. Moreover we have specified
that the order parameter range to consider be split into 100 bins/states (via the first argument); but curiously we have chosen to consider
an order parameter range between 1.0E10 and 1.0E10 (second and third arguments). In FED calculations moves which take the system outwith the
specified range are rejected; the system is constrained to reside within the specified range. Since the aim of our preliminary simulations
is to probe what would be a reasonable range to specify for later FED calculations, we do not wish to artificially constrain the range here, which is
why we have set the range to be 'infinite'. In a similar vein, while we choose in line 10 to use the expanded ensemble
method (described later) to learn the bias function, we specify an update time larger than the length of the simulation (namely, 1,000,000,000 moves)
to essentially 'switch off' this functionality for now, ensuring that the simulation is unbiased and thus visits the equilibrium
configurations as would normally be the case. (Of course one does not necessarily have to use the expanded method at this
point. Any method for learning the bias function would be fine to specify with *fed method* here  so long as biasing is essentially
switched off).
The above 'hacks' may all seem a little strange, but it is currently the only way in DL_MONTE to
track the LSMC order parameter while performing a 'conventional' Monte Carlo simulation, i.e. a Monte Carlo simulation without using lattice switch
moves or biasing  as is required in our preliminary simulations.
*Exercise*

The input files correspond to a preliminary simulation for phase 1. Create a directory in which to perform this simulation, and
copy the input files into it. Run the simulation. The simulation should take about 90 seconds to complete.
Now create another directory for the phase2 preliminary simulation, and copy the input files into it. Change *initactive* in
CONTROL so that the simulation is instead performed in phase 2, and run the simulation. Again, the simulation should take
about 90 seconds to complete.
Postprocessing YAMLDATA

Each preliminary simulation will output a file YAMLDATA.000, which contains a header of metadata (e.g. the temperature and pressure
of the system), followed by data output periodically during the simulation (in YAML format).
(The frequency of output to YAMLDATA.000 is controlled by the *yamldata* directive in CONTROL).
Here are the first few lines of a YAMLDATA.000 output by an instance of the above phase1 preliminary simulation:
.. codeblock:: html
:linenos:
temperature: 1.00000000000000E01
pressure: 0.00000000000000E+00
usingbias:
usingpsmc:


timestamp: 216
bias: 0.00000000000000E+00
orderparam: 2.30319435732781E01
psmcactive: 1
energy: 1.80527415326411E+03
enthalpy: 1.80527415326411E+03
energyvdw: 1.80527415326411E+03
volume: 1.96415444935759E+02

timestamp: 432
bias: 0.00000000000000E+00
orderparam: 2.29806797369747E01
psmcactive: 1
energy: 1.80506844807347E+03
enthalpy: 1.80506844807347E+03
energyvdw: 1.80506844807347E+03
volume: 1.96357624881152E+02

timestamp: 648
bias: 0.00000000000000E+00
orderparam: 2.33436638488456E01
psmcactive: 1
energy: 1.80424889431791E+03
enthalpy: 1.80424889431791E+03
energyvdw: 1.80424889431791E+03
volume: 1.96387869892890E+02

The meaning of the data in the file should be transparent, though note that *psmcactive:* corresponds to the phase which
the current configuration belongs to. We are interested in the range of :math:`M` exhibited by the system in each
phase, which corresponds to *orderparam:* in YAMLDATA.000. The script *strip_yaml.sh* in the *script* directory can be
used to postprocedd YAMLDATA.000, extracting data from it and plotting it against simulation time (i.e. the number of moves so far).
The argument to the script determines which piece of data to extract, as well as the name of the file to store the extracted
data in. For example ``strip_yaml.sh energy`` would create a file energy.dat which contains the energy of the system
vs. time.
*Exercise*

Use *strip_yaml.sh* to create files energy.dat and volume.dat for the preliminary simulations for each phase, which contain, respectively, the energy
vs. simulation time and the volume vs. simulation time for the preliminary simulations. Plot the files, and check that the preliminary simulations
are wellequilibrated within the simulation time.
Example plots of energy.dat and volume.dat from both preliminary simulations are given below.
.. figure:: ImagesPSMC/tutorialpsmc_prelimenergy.jpg
:width: 600px
Energy vs. simulation time for preliminary simulations. Note that the energy equilibrates quickly, within the first 50,000
moves. Note also that the energies of the hcp and fcc phases are very similar.
.. figure:: ImagesPSMC/tutorialpsmc_prelimvolume.jpg
:width: 600px
Volume vs. simulation time for preliminary simulations. Note that the volume equilibrates quickly, within the first 50,000
moves. Note also that the volumes of the hcp and fcc phases are very similar.
Use *strip_yaml.sh* to create files orderparam.dat for the preliminary simulations for each phase, which contain :math:`M` vs. simulation
time for the preliminary simulations.
Plot the files and come up with a range of :math:`M` which encompasses the equilibrium configurations for both phases.
An example plot of orderparam.dat for both preliminary simulations is given below.
.. figure:: ImagesPSMC/tutorialpsmc_prelimorderparam.jpg
:width: 600px
Order parameter vs. simulation time for preliminary simulations. Note that the order parameter equilibrates quickly,
within the first 50,000 moves. Note also that the order parameter distinguishes the hcp and fcc phases: the order parameter is
negative for hcp configurations, while it is positive for fcc configurations.
Discussion

From the above figure it is clear that the hcp phase corresponds to an order parameter range of about 11 to 3, and the fcc
phase corresponds to a range of about 2 to 11. Thus the range 11 to 11 encompasses both phases. However we will use the larger
range 15 to 15 be safe. Configurations with order parameters less than 11 or greater than 11 may occur rarely at equilibrium
though not be visible on our plot. However they still may make a significant contribution to the hcpfcc free energy difference.
However we expect that configurations with order parameters outwith the range 15 to 15 will be so rare to be insignificant.
Biasgeneration simulation
==========================
Having determined an appropriate range for the order parameter (15 to 15), we are now in a position to perform the
biasgeneration simulation.
CONTROL

For this simulation the files FIELD, CONFIG, CONFIG.1 and CONFIG.2 are identical
to the preliminary simulations. However the CONTROL file is slightly different:
.. codeblock:: html
:linenos:
PSMC simulation LJ # Comment line
use fed psmc # Use PSMC in 'FED' functionality, followed by block of PSMC flags; begins FED block
switchfreq 216 # Latticeswitch move frequency
initactive 1 # Phase to start in (1 or 2; chooses configuration in CONFIG)
datafreq 43200 # Output frequency to PSDATA* (default = 1000)
psmc done # End of PSMC block
fed method ee 1.0 1.0 43200000 # Use EE method for learning the bias function; last argument is update frequency
fed order psmc 100 15.0 15.0 1 # Bias over order parameter for PSMC in specified range
fed done # End of FED block
finish # End of 'use' block
ranseed # Use a random seed
steps 432000000 # Simulation length
temperature 0.1
pressure 0.0
sample coordinate 4320000 # Frequency to store configurations during simulation
archiveformat dlpoly4 # Configurations during simulation stored in DL_POLY4 format (HISTORY)
yamldata 43200 # Frequency to store energies etc. in yaml format (YAMLDATA)
move atoms 1 216 # Move atoms, followed by frequency
LJ core
move volume ortho log 1 # Move volume, followed by freqency
print 43200000 # Output frequency to OUTPUT.000 (not needed here  set frequency large)
stat 43200000 # Output frequency to PTFILE.000 (not needed here  set frequency large)
start
The important features of this file are as follows. Firstly, *switchfreq* (line 5) is no longer 0; lattice switch moves are enabled.
Note that the frequency of lattice switch moves is high: here we attempt lattice switch moves with the same frequency as atom translation
moves (i.e. with frequency 216). Because of the way LSMC is implemented in DL_MONTE, there is essentially no additional computational cost
associated with attempting lattice switch moves. (Without going into detail, the computational cost of calculating the energy of the trial
configuration generated by a lattice switch move is, in fact, absorbed into the calculation of the LSMC order parameter).
Hence one may as well attempt lattice switch moves frequently, which ultimately improves the rate of successful transitions between phases.
Secondly, we have set the order parameter range to be from 15 to 15 (line 13). Moreover we have divided the order parameter range into
100 'states'  the grid on which the bias function is defined.
Thirdly, we have now 'enabled' biasing (line 12), via the expanded ensemble method. This is the simplest method for learning the bias function.
In the expanded ensemble method the bias is updated after every block of :math:`n` moves. During each block we track how many the times the
system visits each order parameter state. This information is then used to generate the bias function for the next block using the following
equation:
.. math::
\eta_i^{(k+1)}= \eta_i^{(k)}+\ln(h^{(k)}_i)
where :math:`\eta_i^{(k)}` is the bias for state :math:`i` during block :math:`k`, and :math:`h^{(k)}_i` is the number
of visits to state :math:`i` during block :math:`k`. This equation increases the bias for states which were visited rarely during
the last block, ensuring that they are visited more often during the next iteration. The end result is that after enough blocks all states are sampled
equally. The forthcoming exercise will illustrate this. Here we have chosen :math:`n` to be 43,200,000, which corresponds to every 100,000
'sweeps' of the system (where sweep here is comprised with, on average, 216 atom translation moves, 216 lattice switch moves, and 1 volume
move). This is the origin of the third argument to *fed method ee* on line 12. Note that it is necessary for :math:`n` to be large in the
expanded ensemble method  the block must be large enough that a representative histogram :math:`h^{(k)}_i` is obtained for the current bias function
:math:`\eta_i^{(k)}`. (Do not worry about the first two arguments to *fed method ee*, which are both '1.0'. By modifying these you can alter the
form of the above equation. However this is unimportant here).
Fourthly, we have set the length of the simulation to be 432,000,000 moves; it can take a long time for an accurate bias function to be obtained.
Note that this corresponds to 10 updates of the bias function.
Finally, though this is a minor point, we have used the *sample coordinate*, *print* and *stat* directives so that we only occasionally output
data to the files ARCHIVE.000, OUTPUT.000 and PTFILE.000  in order to keep the size of these files, which we do not use here, small.
The directive *datafreq* in the PSMC block has the same purpose for the LSMCspecific output file PSDATA.000, which we also do not use here.
*Exercise*

Create a new directory in which to perform the biasgeneration simulation. Copy the FIELD, CONFIG, CONFIG.1 and CONFIG.2 input files for one of the
preliminary simulations into this directory, and create a new CONTROL file with the content given above. Run the simulation, noting that it should
take about 56 hours to complete.
The first and third columns in the files FEDDAT.000_001, FEDDAT.000_002, ..., FEDDAT.000_010 contain histograms over :math:`M` (i.e. how often
each order parameter state was visited) for blocks 1, 2, ... 10 during the simulation.
Plot the first and third columns in the FEDDAT file (e.g. using the command
.. codeblock:: html
plot "FEDDAT.000_001" using 1:3,"FEDDAT.000_002" using 1:3",...
in gnuplot) to inspect how the range of :math:`M` explored during each block changes as the block number (i.e. the index on the FEDDAT file name)
increases.
You should see that the explored range increases with block number, until eventually the whole range is explored during a single block, as in the example
plot below:
.. figure:: ImagesPSMC/tutorialpsmc_biasgen_hist.jpg
:width: 600px
Orderparameter histograms for various blocks during a biasgeneration simulation using the expanded
ensemble method. Note that the histograms become increasingly flat as the block number is increased, and for
blocks late in the simulation the entire orderparameter range is explored.
As you did previously during the preliminary simulations, use the *strip_yaml.sh* script to create a file orderparam.dat containing the
order parameter vs. simulation time for the biasgeneration simulation, and plot this. You should see that in the early stages of the simulation
the order parameter is confined to small ranges of :math:`M`, but eventually the simulation sweeps quickly across the entire range of :math:`M` 
in a manner which reflects your previous plot. The orderparam.dat corresponding to the above plot is given below:
.. figure:: ImagesPSMC/tutorialpsmc_biasgen_orderparam.jpg
:width: 600px
Order parameter vs. simulation time for a biasgeneration simulation using the expanded ensemble method.
The first and second columns in the files FEDDAT.000_001, FEDDAT.000_002, ..., FEDDAT.000_010 contain the bias functions calculated at the end of blocks
1, 2, ..., 10.
Plot the first two columns in the FEDDAT files to inspect how the bias function evolves throughout the simulation. You should see that it converges
towards a doublewell shape, as in the example plot below:
.. figure:: ImagesPSMC/tutorialpsmc_biasgen_bias.jpg
:width: 600px
Evolution of the bias function during a biasgeneration simulation using the expanded ensemble method.
Discussion

A wellconverged bias function tells us a lot about the stability of the two phases in question; we need only actually proceed to the production simulation
if we want quantitative results. It turns out that if the bias function is 'perfect', i.e. it leads to perfectly uniform sampling over the order parameter
range, then the bias function is equivalent to the free energy profile over that order parameter.
With this in mind, consider the final bias function obtained from the above simulation.
Clearly it has two minima, one at :math:`M\approx 5` and one at :math:`M\approx 5`.
These are free energy minima corresponding to the equilibrium configurations of each phase: the :math:`M\approx 5` corresponds to the hcp equilibrium
configurations while the :math:`M\approx 5` corresponds to the fcc equilibrium configurations. Between these minima is a maxima. This is the free energy
barrier separating the two phases in the LSMC calcualtion. Note that the fcc minimum is lower than the hcp minimum, suggesting that the fcc phase is prefered
at this temperature and pressure.
However, the above analysis supposes that the bias function obtained from the biasgeneration simulation is perfect. It is not.
Above we noted that the bias function is well converged. However it is not perfectly converged.
Nevertheless at this point we can be fairly confident that the fcc phase is stable over hcp at this temperature and pressure.
There is still the question though of what the precise value of the free energy difference is. We address this with the production simulation.
Biasgeneration simulation (transition matrix)
==============================================
*This section is somewhat of an aside. Feel free to skip this section and proceed to the production simulation, perhaps returning to this section
afterwards if you wish.*
In 'real world' LSMC calculations the biasgeneration simulation is the most technically difficult part. The bias function must be 'good enough' that it
leads to both phases being explored in a reasonable timescale. However it becomes increasingly difficult to generate such a bias function as the system
size increases.
This is due to the fact that the height of the free energy barrier separating the two phases at :math:`M\approx 0`
increases with system size, making it increasingly difficult for the system to explore the entire range of :math:`M`  even with biasing.
It turns out that the method we used above, the expanded ensemble method, is
only tractable for small systems. In this section we repeat the biasgeneration simulation using a superior method, the
*transition matrix method*, and contrast its performance against the expanded ensemble method.
The transition matrix method learns the bias function through tracking the proposed transitions between different order parameter states and their
associated probabilities of being accepted. Details of the method can be found elsewhere (see the links at the end of this tutorial).
As well as being more efficient (at least it is in this simulation) than the expanded ensemble method, the transition matrix method has the further benefit that
it can can exploit parallelisation to reduce the computational time required to generate a good bias function. Multiple simulations can be
run in parallel using the transition matrix method, and their results can be pooled together at their completion
to generate a bias function. This is not possible with the expanded ensemble method. The simulations could even be assigned distinct regions of
the considered order parameter range, reducing computational time further.
CONTROL

The CONTROL file for our biasgeneration simulation exploiting the transition matrix method is as follows (all other input files are unchanged
compared to the previous simulations):
.. codeblock:: html
:linenos:
PSMC simulation LJ # Comment line
use fed psmc # Use PSMC in 'FED' functionality, followed by block of PSMC flags; begins FED block
switchfreq 216 # Latticeswitch move frequency
initactive 1 # Phase to start in (1 or 2; chooses configuration in CONFIG)
datafreq 43200 # Output frequency to PSDATA* (default = 1000)
psmc done # End of PSMC block
fed method tm 4320000 4320 # Use TM method for learning the bias function; output and update frequencies
fed order psmc 100 15.0 15.0 1 # Bias over order parameter for PSMC in specified range
fed done # End of FED block
finish # End of 'use' block
ranseed # Use a random seed
steps 43200000 # Simulation length
temperature 0.1
pressure 0.0
sample coordinate 4320000 # Frequency to store configurations during simulation
archiveformat dlpoly4 # Configurations during simulation stored in DL_POLY4 format (HISTORY)
yamldata 4320 # Frequency to store energies etc. in yaml format (YAMLDATA)
move atoms 1 216 # Move atoms, followed by frequency
LJ core
move volume ortho log 1 # Move volume, followed by freqency
print 43200000 # Output frequency to OUTPUT.000 (not needed here  set frequency large)
stat 43200000 # Output frequency to PTFILE.000 (not needed here  set frequency large)
start
Note that we have enabled the transition matrix method through the *fed method tm* directive on line 12. The first argument to this directive is how
often the bias function is output to a FEDDAT file. Here we choose this to be every 4,320,000 moves. The second argument is how often the bias function
is updated using the information regarding the transitions between order parameter states gathered so far in the simulation.
It is beneficial to perform such updates frequenty, but it is unnecessary (and would result in an unnecessary computational expense)
to do so every move. Here we set the updates to occur every 4320 moves, which corresponds to every 10
'sweeps' of the system (where a sweep here is comprised with, on average, 216 atom translation moves, 216 lattice switch moves, and 1 volume move).
Note also that this simulation is far shorter than the previous biasgeneration simulation: this simulation is of length 43,200,000 moves, while the
previous one using the expanded ensemble method was of length 432,000,000 moves  a factor of 10 longer. We use such a short simulation here in
anticipation of the fact that the transition matrix method is far more efficient than the expanded ensemble method, and hence will yield a good
bias function in far less moves. As we shall see in a moment, 43,200,000 moves is indeed sufficient for the transition method to yield a good bias function,
while this was far from the case using the expanded ensemble.
*Exercise*

Create a new directory in which to perform the biasgeneration simulation using the transition matrix method. Copy the FIELD, CONFIG, CONFIG.1 and CONFIG.2
input files for one of the preliminary simulations into this directory, and create a new CONTROL file with the content given above.
Run the simulation, noting that it should take about 30 minutes to complete.
As you did previously for the biasgeneration simulation using the expanded ensemble method, plot the FEDDAT files and inspect the bias function throughout
this biasgeneration simulation. You should see that the bias function is well converged on the timescale of the simulation, as in the example plot below:
.. figure:: ImagesPSMC/tutorialpsmc_biasgentm_bias.jpg
:width: 600px
Evolution of the bias function during a biasgeneration simulation using the transition matrix method.
Discussion

Note that here the transition matrix method converged to the correct bias function an order of magnitude faster than the expanded ensemble method.
This serves to illustrate that significant speedups can be realised by using a more efficient algorithm.
In fact, the transition matrix method can be exploited to realise even greater speedups than demonstrated here, by exploiting parallelisation. The
relevant functionality exists in DL_MONTE. However it is beyond the scope of the current tutorial to demonstrate it.
Production simulation
=====================
We now turn to the production simulation, using the bias function obtained from the biasgeneration simulation. Recall the aim of this simulation:
to, with the aid of lattice switch moves and biasing, sample from both phases, and measure :math:`p_1/p_2` (where :math:`p_1` and :math:`p_2`
are the probabilities of the system being in phases 1 and 2 respectively), which can then be used to obtain the free energy difference
:math:`\Delta G` via the equation
.. math::
\Delta G \equiv G_1  G_2 = k_BT\log(p_1/p_2)
There is, of course, the question of how we obtain :math:`p_1/p_2` from the production simulation. After all, the production simulation uses biasing:
how often phases 1 and 2 are visited during the production simulation does not reflect how often phases 1 and 2 are visited in the 'true' system (the
quantities :math:`p_1` and math:`p_2` pertain to the true system). This is because we 'force' the simulation to visit both phases approximately
equally using biasing. How often phases 1 and 2 are visited in the simulation in fact reflects the probabilities of the system being in either phase
in the *biased* system. What we must do to obtain :math:`p_1` and :math:`p_2` (and hence :math:`p_1/p_2`) is to take the data from the production
simulation, and *reweight* it so that it pertains to the unbiased system. We will address this later, but first we must run the production simulation...
Input files

The input files for the production simulation are identical as for the biasgeneration, with the following exceptions.
Firstly, we specify the bias function from the outset via an input file FEDDAT.000  as opposed to learning the bias function during the simulation. The
last FEDDAT file output from the biasgeneration simulation, FEDDAT.000_010, contains the bias function we wish to use in the production simulation.
Conveniently, this file can be used, without modification, as the input file FEDDAT.000 for the production simulation in order to specify the bias function.
However one must instruct DL_MONTE that the bias function is to be read from FEDDAT.000, as opposed to being initialised as zero for all order parameter
states as was the case in previous simulations. This is done by changing
.. codeblock:: html
fed order psmc 100 15.0 15.0 1
in CONTROL to
.. codeblock:: html
fed order psmc 100 15.0 15.0 1
where the '1' signifies that the bias function is to be read from FEDDAT.000.
The other change which we must make to the input files, which is also to CONTROL, is that we must ensure that the bias function is not updated during
the simulation, since we wish the bias function to remain fixed throughout. Recall that this was accomplished for the preliminary simulations by using the
directive
.. codeblock:: html
fed method ee 1.0 1.0 1000000000
This directive is also suitable here.
*Exercise*

Create a new directory in which to perform the production simulation. Copy the FIELD, CONFIG, CONFIG.1, CONFIG.2 and CONTROL input files for the
biasgeneration simulation into this directory. Copy the FEDDAT.000_010 output file from the biasgeneration simulation to this directory,
and rename it FEDDAT.000. Modify the CONTROL file as described above so that it corresponds to a production simulation: this entails
modifying the two lines containing the *fed order* and *fed method* directives as described above.
Run the simulation, noting that it should take about 56 hours to complete.
Once the simulation is complete, there will be a file FEDDAT.000_001 which contains the bias function (first and second columns) and a histogram over
:math:`M` for the whole simulation (first and third columns). Note that the bias function will be identical to that in
FEDDAT.000, since we instructed DL_MONTE not to update it during the simulation. As you did earlier for the FEDDAT files output by the biasgeneration
simulation, plot the first and third columns in FEDDAT.000_001 to see the aforementioned histogram. You
should see that the entire order parameter range was sampled approximately uniformly, as was desired  as in the example plot below:
.. figure:: ImagesPSMC/tutorialpsmc_prod_hist.jpg
:width: 600px
Orderparameter histogram for the production simulation.
As you did previously during the preliminary and biasgeneration simulations, use the *strip_yaml.sh* script to create a file orderparam.dat
containing the order parameter vs. simulation time, and plot this. You should see that the simulation efficiently explores the entire order
parameter range, as in the example plot below:
.. figure:: ImagesPSMC/tutorialpsmc_prod_orderparam.jpg
:width: 600px
Order parameter vs. simulation time for a production simulation.
Reweighting to remove the bias

Above you plotted the order parameter histogram for the production simulation. This information was contained in the first and third columns of the file
FEDDAT.000_001. Let :math:`h_i` denote how often order parameter state :math:`i` was visited during the simulation.
Conveniently, the bias :math:`\eta_i` for each order parameter state :math:`i` is
contained in the second column of FEDDAT.000_001. Now, the effect of biasing is to add an energy :math:`\eta_i/\beta` to all configurations belonging to
order parameter state :math:`i`. This amplifies the sampling of state :math:`i` by a factor :math:`\exp(\eta_i)`: if :math:`i` were visited
with frequency :math:`f` during an unbiased simulation, it would be visited with frequency :math:`f\exp(\eta_i)` in a biased simulation. Conversely
if :math:`i` were visited with frequency :math:`f'` in a biased simulation, it would be visited with frequency :math:`f'\exp(\eta_i)` in
an (ergodic) unbiased simulation. Thus we can obtain the orderparameter histogram corresponding to an unbiased simulation by multiplying
the frequency corresponding to each state in the biased simulation by :math:`\exp(\eta_i)`.
*Exercise*

The script *unbiasedhist.sh* performs the aforementioned reweighting, using the data in FEDDAT.000_001 file in the current directory. It
creates a file unbiasedhist.dat which contains the corresponding unbiased histogram. Apply this script to the output of your production simulation,
and plot the resulting file unbiasedhist.dat. An example plot is given below:
.. figure:: ImagesPSMC/tutorialpsmc_prod_unbiasedhist.jpg
:width: 600px
*Unbiased* orderparameter histogram created by reweighting the biased orderparameter histogram obtained from the production simulation.
Note that there are two peaks in the unbiased histogram. These correspond to the equilibrium configurations for hcp (left) and fcc (right) phases. Note
also that the fcc peak is larger than the hcp peak, signifying that the fcc phase is stable over hcp at this temperature and pressure. Thus what we
saw previously earlier for the bias function while analysing our data from the biasgeneration simulation is borne out in the data from the production
simulation.
Recall that negative order parameters correspond to phase 1 while positive order parameters correspond to phase 2. With this in mind, the unbiased
histogram can be used to obtain :math:`p_1/p_2` as follows: sum all elements of the histogram corersponding to phase 1, do the same for elements corresponding
to phase 2, and then divide the former by the latter. The script *unbiasedhist_probratio.sh* performs this task, acting on the unbiasedhist.dat file in the
current directory. It outputs a single value, :math:`p_1/p_2`.
Apply *unbiasedhist_probratio.sh* to extract :math:`p_1/p_2` for your production simulation. Substitute this quantity into the expression for
:math:`\Delta G` given at the beginning of this section in order to obtain a value for :math:`\Delta G`, remembering that our simulations used
:math:`T=0.1\epsilon/k_B`. Divide this value of :math:`\Delta G` by the number of particles in the system, :math:`N`, to get the free energy difference
per particle :math:`\Delta G/N` (in units of :math:`\epsilon`).
The 'true' value of :math:`\Delta G/N` for this system, obtained using long DL_MONTE simulations and verified against another LSMC code, is
:math:`\Delta G/N=0.001283(7)\epsilon`, where the '7' in parenthesis is the uncertainty on the final digit (the standard error of the mean  see below).
Is your value of :math:`\Delta G/N` in agreement with the true value?
Calculating uncertainties

Above we evaluated :math:`\Delta G` by deducing an unbiased orderparameter histogram for the simulation, integrating over the regions
of the histogram corresponding to phases 1 and 2 to obtain :math:`p_1/p_2`, and then substituting this :math:`p_1/p_2` into the equation relating
:math:`p_1/p_2` to :math:`\Delta G`. This procedure certainly yields a value for :math:`\Delta G`. However there is the question of what the
*uncertainty* in the :math:`\Delta G` obtained from our production simulation is. The above method for obtaining :math:`\Delta G` does not provide this.
Moreover the histogram used in the analysis above takes into account the configurations visited at the very start of the simulation  during the
equilibration time. These configurations should not be considered when evaluating :math:`\Delta G` (or indeed any other quantity). Here we
demonstrate how to calculate :math:`\Delta G` and its associated uncertainty more rigourously using the data output by the production simulation.
We first briefly revise how to combine multiple measurements of a physical quantity :math:`X` into a single 'quoted value' and an associated uncertainty.
Let us assume that we have measured :math:`X` :math:`n` times. Let :math:`X_1,X_2,\dotsc,X_n` denote these :math:`n` measurements.
The standard practice is to calculate the mean of all the measurements,
.. math::
\bar{X}=\frac{1}{n}\sum_{i=1}^nX_i,
and then use this as the 'quoted value' of :math:`X`. Regarding the uncertainty in :math:`\bar{X}`, one commonly uses the *standard error of the mean*:
.. math::
\delta X=\text{stdev}(X)/\sqrt{n}
as the uncertainty in :math:`\bar{X}`, where :math:`\text{stdev}(X)` is the standard deviation of the observations :math:`X_1,X_2,\dotsc,X_n`.
However it is also common to use *twice* the standard error of the mean as the uncertainty. (This corresponds to a 95% confidence interval, while
the standard error of the mean corresponds to a 68% confidence interval).
Due to the use of random numbers in Monte Carlo, diffferent simulations (using different randomnumbergenerator seeds) will yield different
values for a physical quantity :math:`X`. We can treat the results of different simulations as different 'measurements' in the sense above,
and combine the results from these simulations to obtain a quoted value and an associated uncertainty as described above.
In fact, one does not need to resort to combining the results from separate simulations. Instead, one can split the data obtained from a single
simulation into blocks, calculate :math:`X` for each of these blocks, and use these :math:`X` as the measurements from which the quoted value
and uncertainty are calculated.
We apply this approach below to the data output by the production simulation to calculate a quoted value and associated uncertainty
for :math:`\Delta G`.
*Exercise*

The script *probratio.sh* processes the data in YAMLDATA.000 to calculate :math:`p_1/p_2`, but only the data in YAMLDATA.000 which corresponds to
a specified 'time window'. The script takes two arguments, :math:`t_{\text{min}}` and :math:`t_{\text{max}}`, which correspond to the lower and
upper bounds of the window  where the units of time here are the number of elapsed Monte Carlo moves (which corresponds to the *timestamp:*
data in YAMLDATA.000). An example usage is as follows, where the script is called in the directory containing the YAMLDATA.000 file for the production
simulation (the ``$`` is the shell prompt):
.. codeblock:: html
$ probratio.sh 4320000 432000000
0.0610501
In this example :math:`t_{\text{min}}=4320000` and :math:`t_{\text{max}}=432000000`, and 0.0610501 is the value of :math:`p_1/p_2`
corresponding to the this time window.
Recall that the length of our production simulation was 432,000,000 moves. This can be split into 4 blocks of 108,000,000 moves each,
in which case *probratio.sh* can be used to obtain :math:`p_1/p_2` obtained for each block via
.. codeblock:: html
probratio.sh 1 108000000
for the 1st block,
.. codeblock:: html
probratio.sh 108000001 216000000
for the 2nd block,
.. codeblock:: html
probratio.sh 216000001 324000000
for the 3rd block, and
.. codeblock:: html
probratio.sh 324000001 432000000
for the 4th block. The :math:`p_1/p_2` obtained for each block in this manner can then be used to obtain a :math:`\Delta G` for each
block as described earlier, after which the mean and standard error of the mean of all blocks can be calculated.
Do this for your production simulation, and compare your result to the 'true' value given above.
Do the same again, but this time assume an equilibration time of 4,320,000 moves, i.e. ignore the first 4,320,000 moves, partitioning
the remaining time into four blocks.
One concern is how large each block should be. The answer is that the 'block size', i.e. the size of the time window each block coversm should be
larger than the correlation time for the physical quantity :math:`X` one is interested in. This ensures that the blocks, and hence also their
associated values of :math:`X`, are uncorrelated. For the purpose of evaluating :math:`\Delta G` the block size should be larger than the time
it takes for the system to explore both phases. One can check this by plotting the order parameter against simulation time, and inspecting the
timescale of the fluctuations in the order parameter against the block size: the timescale of the fluctuations should be far smaller than the block
size. Verify that this is indeed the case here. (Recall that you have already plotted the order parameter against simulation time for the
production simulation earlier).
As one final check, repeat the above for various block sizes. The essential result should be independent of block size  so long as one does
not make the block size so small that both phases are not explored in any block, or that the block size is so small that time resolution of the
data in YAMLDATA.000 is too coarse to yield meaningful results.
Further information
===================
 Section 2 of the paper at https://arxiv.org/abs/1609.04329 (which describes an old LSMC code) provides a more thorough discussion of LSMC relevant
to this tutorial, including the expanded ensemble method (called the 'visited states method' there), the transition matrix method, and a
description of the LSMC order parameter :math:`M` relevant to the simulations performed here
 The DL_MONTE reference manual provides more detailed information regarding LSMC functionality in DL_MONTE